#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<ctime>

using namespace std;

#define INTERATIONS 100000000


long double borwein(int N){

	long double a, b, p;

	a = sqrt(2);
	b = 0;
	p = 2 + sqrt(2);

	for(int i = 0; i < N; i++){
		b = (sqrt(a) * (1 + b)) / (a + b);
		a = ((sqrt(a)) / 2) + (1 / (2 * sqrt(a)));
		p = p * b * ((1 + a) / (1 + b));
	}
	
	
	return p;
}

long double gauss(int N){

	unsigned int p;
	long double a0,a1, b, t, pi;
	
	a0 = 1;
	b = 1/sqrt(2);
	t = 0.25;
	p = 1;
	
	for(int i = 0; i < N; i++){
		a1 = (a0+b)/2;
		b = sqrt(a0 * b);
		t = t - p * pow((a0 - a1),2);
		p = 2*p;
		pi = pow((a1 + b),2)/(4 * t);
		a0 = a1;
	}
	
	return pi;
}

double mean(double *vet){
	int size = sizeof(*vet)/sizeof(double);
	//cout << size << endl;
	double sum = 0;
	
	for(int i = 0; i < size; i++){
		sum += vet[i]/size;
	}
	
	return sum;
	
}

double stddev(double *vet,double mean){
	double sum = 0;
	int size = sizeof(vet)/sizeof(double);
	
	for(int i = 0; i < size; i++){
		sum += pow((vet[i] - mean),2) / (size - 1);
	}
	
	
	return sqrt(sum);
}

long double monte_carlo(){

	float S,
	E,
	r,
	sig,
	T;
	int in = 0;	// points in

	double pi,
	dist,	// distance
	x,	// x value
	y;	// y value

	srand(time(NULL));

	for(int i = 0; i < N; i++){
		x = (double) rand()/RAND_MAX;
		y = (double) rand()/RAND_MAX;
		dist = sqrt(pow(x,2) + pow(y,2));

		if(dist <= 1)	in++;
	}
	pi =(double) 4 * in/N;
	
	return pi;
}



int main(int argc,char**argv){
	
	if(argc < 2){
		cout << "usage ./prog <option>" << endl;
		cout << "option = 1 - borwein" << endl;
		cout << "option = 2 - gauss" << endl;
		cout << "option = 3 - monte-carlo" << endl;
	}
	
	long double pi = 0;
	
	
	int option = atoi(argv[1]);
	
	if(option == 1){
		pi = borwein(INTERATIONS);
	} else if(option == 2){
		pi = gauss(INTERATIONS);
	} else if(option == 3){
		pi = monte_carlo();
	} else{
		cout << "invalid option: choose '1', '2' or '3'" << endl;
		return -1;
	}
	
	cout.precision(7);
	cout << "pi:" << pi << endl;
	
	
	return 0;
}
